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Abstract 

We have already reported the first result on the all-particle spectrum around the 
knee region based on data from 2000 November to 2001 October observed by the 
Tibet-Ill air-shower array. In this paper, we present an updated result using data set 
collected in the period from 2000 November through 2004 October in a wide range 
over 3 decades between 10 14 eV and 10 eV, in which the position of the knee is 
clearly seen at around 4 PeV. The spectral index is -2.68 ± 0.02(stat.) below IPeV, 
while it is -3.12 ± 0.01(stat.) above 4 PeV in the case of QGSJET+HD model, and 
various systematic errors are under study now. 
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1 Introduction 



The energy spectrum of primary cosmic rays is well described by a power law 
over a wide energy range, while its slope becomes steeper in the energy range 
between 10 15 and 10 16 eV, which is called "knee". It has been discussed that 
the knee is closely related with the origin, the acceleration and the propagation 
of high-energy cosmic rays in the Galaxy. One of the plausible understanding 
may be that almost all cosmic rays below the knee are accelerated at supernova 
remnants (SNRs), since the maximum energy gained by shock acceleration at 
SNRs is of the order of 10 14 eV per unit charge [Lagage and Cesarsky. 1983], 
the cosmic-ray spectrum is expected to become steeper at energies around and 
beyond the knee. Another possibility is that the break of the spectrum around 
the knee represents the energy at which cosmic rays can escape more freely 
from the trapping zone in the Galactic disk [Peters 1961]. 

We have already reported the first result on the all-particle spectrum around 
the knee region based on data from 2000 November to 2001 October observed 
by the Tibet-Ill air-shower array [Amenomori et al., 2003]. In this paper, we 
present an updated result using data set collected in the period from 2000 
November through 2004 October. We also examine the simulation code COR- 
SIKA with interaction models of QGSJETOlc. Based on these calculations, 
we obtained the cosmic ray energy spectrum in a wide range over 3 decades 
between 10 14 eV and 10 17 eV. 



2 Tibet-Ill air shower array and detector response 

The Tibet-Ill air-shower array, consisting of 533 scintillation detectors of each 
0.5 m 2 with the area of 22,050 m 2 , has been successfully operating since 1999 
with energy threshold of a few TeV. A 0.5 cm-thick lead plate is put on the top 
of each detector to increase the detection sensitivity of a detector by converting 
secondary gamma rays into electron-positron pairs in an air shower. In the fall 
of 2000, this array was further improved for UHE cosmic-ray study by adding 
wide dynamic range PMTs to all 221 sets of detectors which are placed on a 
lattice of 15 m spacing in the detector covering area of 36,900 m 2 . This PMT 
equipped in each detector can measure the number of particles beyond 4,000, 
so that the array can observe UHE cosmic rays with the energy exceeding 10 17 
eV with a good accuracy. 

The Tibet-Ill air-shower array is able to measure the shower size and the 
arrival direction of each air shower. The air shower direction can be estimated 
with an inaccuracy smaller than 0.2° at energies above 10 14 eV. The shower size 
N e for each event is estimated by fitting the lateral particle density distribution 
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with the modified NKG structure function (see section 3.2). The primary 
energy of each air shower event is then estimated from the air shower size. 
The relation between the shower size and the primary energy is calculated 
based on the MC simulation as discussed later. In general, the "shower size" 
means the number of electrons and positrons in each air shower event. It cannot 
be, however, directly obtained by usual shower array using scintillators since 
these electromagnetic components cannot be discriminated from other charged 
particles such as muons and hadrons. 

In our experiment, the number of particles detected by each scintillation de- 
tector is defined as the PMT output (charge) divided by that of the single 
peak, which is determined by a probe calibration using cosmic rays. For this 
purpose, a small scintillator of 25 cm x 25 cm x 3.5 cm thick with a PMT 
(H1949) is put on the top of the each detector during the maintenance pe- 
riod. This is called probe detector and is used for making the trigger of the 
each Tibet-Ill detector. These events triggered by the probe detector was also 
examined by a MC simulation. In this simulation, the primary particles were 
sampled in the energy range above the geomagnetic cutoff energy at Yangba- 
jing, and all secondary particles which pass through the probe detector and 
the Tibet-Ill detector were selected for the analysis. Since the value of PMT 
output is proportional to the energy loss of the particles passing through the 
scintillator, the peak position of the energy loss distribution is defined as the 
"energy deposit by a singly charged particle". This value corresponds to the 
observed single peak of the probe calibration. The peak energy was calculated 
as 6.11 MeV for the Tibet-Ill detector configuration. We confirmed that the 
shape of the energy loss distribution shows a good agreement with the charge 
distribution of the experimental data as shown in Fig. 1. It should be noted 
that the dependences of the energy loss on particle type and its energy are 
adequately taken into account in this calibration. 

For air-shower MC calculation at high energies, we treat the number of shower 
particles as the calculated energy deposit divided by 6.11 MeV. Thus, all de- 
tector responses including muons and the materialization of photons inside 
the detector are taken into account. The shower size of each event was esti- 
mated using a modified NKG lateral distribution function which is tuned to 
reproduce the above defined lateral distribution using the MC events under 
our detector configurations. 



3 Air shower size and primary particle energy 

An extensive Monte Carlo simulation (MC) using CORSIKA code (ver. 6.204) 
including QGSJETOlc [Heck et al., 1998] hadronic interaction model was made 
to obtain the primary cosmic ray spectrum using the Tibet-Ill air shower data. 
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Since the Tibet hybrid experiment of the air shower array and the burst de- 
tector array to measure the energy spectrum of the light components (proton 
and helium) strongly suggests that the knee region is dominated by heavy 
components [Amenomori et al., 2006], in this paper we use a heavy dominant 
(HD) composition model [Amenomori et al., 2000] in the MC for comparison 
with experimental data. All secondary particles are traced until their energies 
become 1 MeV in the atmosphere. Simulated air-shower events were input to 
the detector with the same detector configuration as the Tibet-Ill array with 
use of Epics code (ver. 8.64) [Kasahara] to calculate the energy deposit of 
these shower particles. 

3.1 Event selection and collecting area 

Following criteria are applied to select the events for the analysis. 

1) The zenith angle (9) of each air-shower event should be smaller than 25°, 
or sec 9 < 1.1 to exclude the zenith angle dependence of the air-shower devel- 
opment. 

2) More than 10 detectors should detect a signal of more than five particles 
per detector. 

3) The central positions weighted by the 8th power of the number of particles 
at each detector should be inside the innermost 135 m x 135 m area. This 
area is chosen with use of MC events so that the following two cases are just 
canceling out each other. Namely, the number of events originally inside of 
this area falling outside of this area after event reconstruction equals to the 
number of events in the opposite case. 

It is confirmed by simulations that the air showers induced by primary particles 
with E > 100 TeV and sec 9 < 1.1 can be fully detected without any bias 
under above mentioned criteria. 

The total effective area S x f2 is calculated to be 10410.1 m 2 -sr for all primary 
particles with E > 100 TeV. For the operation period from 2000 November 
through 2004 October, the effective live time T is 2.21 years. The total number 
of air showers selected under the above conditions is 4.1 x 10 7 events. 

3.2 Lateral distribution of shower particles 

The determination of the lateral distribution function of an air shower is sub- 
stantially important in this experiment, since the total number of shower par- 
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tides is obtained by fitting the structure function to the experimental data. 
Since our air-shower array is also used to study the TeV gamma-ray point 
sources, lead plate of 5 mm thick is placed above each scintillation counter to 
increase the sensitivity, the detector response should be carefully examined by 
MC. Using the Monte Carlo data obtained under the same conditions as the 
experiment, we find that the following modified NKG function can fit well to 
the lateral distribution of shower particles under the lead plate: 



C(s) [ r m ' } {+ r m ' ) ,Tm 



C{s) = 2irB(a{s) + 2, -b(s) - a(s) - 2) (2) 



where r m ' = 30 m, and the variable s corresponds to the age parameter, N e 
the total number of shower particles and B denotes the beta function. The 
functions a(s) and b(s) are determined as follows. In CORSIKA simulation, 
the shower age parameter s is calculated at observation level by fitting to a 
function for the one dimensional shower development. It may be possible to as- 
sume that air showers with the same shower age s are in the almost same stage 
of air shower development in the atmosphere, i.e. they show the almost same 
lateral distribution for shower particles irrespective of their primary energies. 
The lateral distribution of the particle density obtained by the simulation with 
carpet array configuration is normalized by the total number of particles which 
is derived from the total energy deposit in infinitely wide scintillator. These 
events are then classified according to the stage of air shower development 
using the age parameter and they are averaged over the classified events. The 
fitting of the equation (1) to the averaged MC data is made to obtain the 
numerical values a and b. Thus, we can obtain the behavior of a and b as a 
function of s as shown in Fig. 2 where original definitions of a(s) and b(s) in 
NKG function are shown by the dotted lines. Although our result shows dif- 
ferent dependences of a and b on s, it is confirmed that the lateral distribution 
of the shower particles is better reproduced by our formula. This expression 
is valid in the range of s = 0.6 ~ 1.6, sec 9 < 1.1 and r = 5 ~ 3000 m. 

Based on the Monte Carlo simulation, the correlation between the true shower 
size and the estimated shower size is demonstrated in Fig. 3(a) and Fig. 3(b). 
As seen in these two figures, a good correlation between the true shower size 
and the estimated shower size is obtained and the shower size resolution is 
estimated to be 7% above N e > 10 5 with sec9 < 1.1. 
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3.3 Relation between the shower size and the primary energy 



In Fig. 4, we show the correlation between the shower size N e under the 
lead plate at the Tibet observation level and the primary energy E for the 
QGSJET+HD model. The conversion function from the shower size N e to the 
primary energy E can be expressed by the following equation for sec 9 < 1.1, 

E = ax(-^-) [TeV], (3) 



where a=1.88, b=0.92. The energy resolution is estimated by our simulation 
as 36% and 17% at energies around 200 TeV and 2000 TeV, respectively. It is 
seen that the energy estimation of the primary particle from the air shower size 
is well made with a good accuracy in the wide energy range including the knee 
region. We further examine the dependence of the result on the interaction 
models, primary models, and other possible systematic biases in the very near 
future. 



4 Results and Discussions 



We obtained the all-particle energy spectrum between 1 x 10 14 eV and 1 
x 10 17 eV based on the QGSJET+HD model as shown in Fig. 5. The red 
circle represents this work, and the blue inverse-triangle is our old result 
[Amenomori et al., 2003], and they are compared with other experiments. Our 
new result shows a good agreement with the previous one in the energy range 
less than 4 PeV, while suggesting a slightly steeper power index above 4 PeV 
and about 10% lower intensity around 10 16 eV. This difference may be due to 
the upgrade of MC calculations and the increase of the observed data. The 
spectral index of this work is -2.68 ± 0.02(stat.) below 1 PeV and -3.12 ± 
0.01(stat.) above 4 PeV in the case of QGSJET+HD model. The evaluation 
of various systematic errors is under examination at present. 

The all-particle spectrum obtained in a wide range over 3 decades between 
10 14 eV and 10 17 eV clearly shows the position of the knee being at energy 
around 4 PeV. In this work, we obtained the result based on the assump- 
tion that the knee energy region is dominated by nuclei heavier than helium 
[Amenomori et al., 2006]. It is noted that the all-particle spectrum depends 
slightly on the primary composition in the energy region below about 5 x 10 14 
eV. A further study of the dependence of the spectrum on the simulation 
codes, interaction models, primary composition models, and other possible 
systematic biases should be carefully done in the very near future. 
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Fig. 1. Distribution of the output charge of PMT by probe calibration of the Ti- 
bet-Ill detector. In order to compare with the simulation on the energy deposit, 
the MC result is adjusted by multiplying a constant to meet with the same peak 
position as the experiment. The fluctuation of the number of photons in scintillation 
light is taken into account with the normal distribution in MC. 
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Fig. 2. The numerical values of a and b are plotted as a function of s, where original 
definitions of a(s) and b(s) in NKG function are shown by the dotted lines, and the 
open circles denote the averaged MC data which are fitted by empirical formulae 
shown by red lines. See the text. 



10 




Fig. 3. (a) The correlation between the true shower size (Netme) and the estimated 
shower size (Ne es t). (b) The shower size resolution is estimated to be 7% above 
iVe > 10 5 with sec/9 < 1.1. 




Fig. 4. Correlations of the air-shower size N e and the primary energy Eq at the 
Tibet observation level for the QGSJET+HD model with secQ < 1.1. The solid 
circles denote the average values, and the solid lines is a fit by equation ( 3 ) 
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Fig. 5. Differential all-particle cosmic-ray flux in a wide range over 3 decades 
between 10 14 eV and 10 17 eV measured by the Tibet-Ill air-shower array 
using the QGSJET+HD model. This work is compared with other experi- 
ments : (JACEE [Asakimori et al., 1998], RUNJOB [Apanasenko et al., 2001], 
PROTON satellite [Grigorov et al., 1971], KASCADE [Antoni et al., 2005], 
BASJE-MAS [Ogio et al., 2004], AKENO [Nagano et al., 1984], CASA-MIA 
[Glasmacher et al., 1999], Tibet-Ill [Amenomori et al., 2003]). 
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